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N-body simulations of coUisionless collapse have offered important clues to the construction of 
realistic stellar dynamical models of elliptical galaxies. Such simulations confirm and quantify 
the qualitative expectation that rapid collapse of a self-gravitating coUisionless system, initially 
cool and significantly far from equilibrium, leads to incomplete relaxation. In this paper we 
revisit the problem, by comparing the detailed properties of a family of distribution functions 
derived from statistical mechanics arguments to those of the products of coUisionless collapse 
found in N-body simulations. 



1 Introduction 



The collapse of a dynamically cold cloud of stars can lead to the formati on o f realistic stellar 
systems, with projected density profiles well represented by the i?^'^ lawl^. The theoretical 
framework for the mechanism of incomplete violent relaxation that governs this scenario of 
structure formation was proposed by Lynden-Bell'^, who argued that fast fluctuations of the 
potential during collapse would lead to the formation of a well-relaxed isotropic core, embedded 
in a radially anisotropic, partially relaxed halo. This general picture served as a physical justi- 
fication for the construction of the so-called /oo models , which indeed recovered the R^''^ law 
and, suitably extended to the case of two-component systems (to account for the coexistence of 
luminous and dark matter), led to a number of interesting applications to the observations'^. 

The application of statistical mechanics to this formation scenariol^ also led to the deriva- 
tion of a separate family of distribution functions, the f^^' models, which was recently shown 
to possess interesting thermodynamic properties in the context of the gravothermal catastro- 
phe"^. The key ingredient for the construction of the f^'^' models is the conjecture that a third 
quantity Q (defined as Q = ji' J'^\E\^'^'^/'^ fSqd^p), in addition to the total mass M and the 



total energy Etot, is approximately conserved during the process of collisionless collapse. This 
quantity is introduced to model the process of incomplete violent relaxation, ensuring a radially 
biased pressure tensor and a 1/r^ density profile in the outer parts of the system. A preliminary 
inspection of the general characteristics of the f^'^' models then convinced us that, with signif- 
icant advantage over the /oo models, they might also serve as a good framework to interpret 
the results of simulations of collisionless collapse not only qualitatively, but also in quantitative 
detail, which is the main point of this paper. 

2 Models 

The extremization of the Boltzmann entropy under the three constraints described in the pre- 
vious section leads to the distribution function f^^> = Aex.p[—aE — d{J'^/\E\^^^y^^], where 
a, A, d, and u are positive real constants; here E {E < 0) and J denote single-star specific 
energy and angular momentum. At fixed value of v, one may think of these constants as pro- 
viding two dimensional scales (for example, M and Q) and one dimensionless parameter, such 
as 7 = ad'^''^ /{AttGA). In the following we will focus on values of u ranging from 3/8 to 1. The 
corresponding models are constructed by solving the Poisson equation for the self-consistent 
mean potential $(r) generated by the density distribution associated with f^'^'. At fixed value 
of z^, the models thus generated make a one-parameter family of equilibria, described by the con- 
centration parameter ^ = — a$(0), the dimensionless depth of the central potential well. The 
projected density profile of concentrated models is well fitted by the R^'"^ law, with n close to 
4, and the models provide reasonable fits to the surface brightness and to the kinematic profiles 
of bright elliptical galaxies'^. 

3 Simulations of collisionless collapse 

In principle, to study the process of collisionless collapse we have two options: either a tree 
codel^or a particle-mesh'^ algorithm. We are interested in the large scale structure of the end- 
products of collisionless collapse, for systems that do not exhibit large deviations from spherical 
symmetry. The natural choice thus appears to be that of a particle-mesh code, based on a 
spherical grid and an e xpansion in spherical harmonics. The code used in the present study is 
thus a new version!^ of the van Albada code. For completeness, we have also run a number 
of comparison simulations with the fast code developed by Dehnenl^, which confirmed that our 
results do not depend on the numerical scheme employed. 

During the process of collisionless collapse, violent relaxation is expected to wipe out much of 
the details that characterize the initial conditions, so that the end-products are expected to have 
properties that depend only or mostly on the initial value of the virial ratio u = {2K/\W\)t=o, 
which sets the relevant collapse factor. In reality, violent relaxation is incomplete. Therefore, 
the final state is that of an approximate dynamical equilibrium characterized by an anisotropic 
distribution function, different from a Maxwellian (which would correspond to thermodynamic 
equilibrium). Because of such incomplete relaxation, the end-products of the simulations do 
conserve some memory of the initial state. 

Earlier investigations' '' compared "clumpy" to "homogeneous" initial conditions, showing 
that clumpy initial conditions lead to end-states with projected density distributions well fitted 
by the R^''^ law. We argue that, if the degree of symmetry in the initial conditions is excessive, 
little room is left for the mechanism of violent relaxation; this is confirmed by the fact that little 
or no mixing is observed in the single-particle angular momentum for homogeneous simulations, 
as reported in Fig. ^ We have thus chosen to focus this paper on clumpy initial conditions, for 
which an efficient mixing in phase space is guaranteed. 
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Figure 1: Single-particle angular momentum scatter (correlation between the initial and final values of J) for 
a typical clumpy simulation (left panel) and for its symmetrized (homogeneous) version (right panel). The 
symmetrization process in the initial state is performed by accepting the radius and the magnitude of the velocity 
of each particle, following the procedure for initial clumpy conditions, but by redistributing uniformly the angular 

variables. 

We performed several runs varying the number of clumps and the initial virial ratio u in 
the range from 0.05 to 0.25. In most cases the clumps are cold, i.e. their kinetic energy is 
all associated with the motion of their centers of mass. The simulations are run for several 
dynamical times beyond the violent collapse phase, ensuring that the system has settled into a 
quasi-equilibrium. The final configurations are quasi-spherical, with shapes that resemble those 
of E2 — E3 galaxies. The final equilibrium half-mass radius r^./ is basically independent of the 
value of u. The central concentration achieved correlates with u, as can be inferred from the 
conservation of maximum density in phase space"^. 



4 Fits and phase-space properties 

In order to study the output of the simulations we compare the density and the anisotropy 
profiles, p{r) and a{r) (here the local anisotropy is defined as a{r) = 2 — {{pg) + {Ps)) / {p1)) ^ o^ 
the end-products with the theoretical profiles of the f^'^' family of models. Smooth simulation 
profiles are obtained by averaging over time in the last few dynamical times of the simulation. 
For the fitting models, the parameter space explored is that of an equally spaced grid in (z/, ^), 
with a subdivision of 1/8 in v, from 3/8 to 1, and of 0.2 in ^, from 0.2 to 14.0; the mass and 
the half-mass radius of the models are fixed by the scales set by the simulations. A minimum 
X^ analysis is performed as described elsewherel^. 

As illustrated in Fig. |21 the density of the simulations is well represented by the best-fit f^'^' 
profile over the entire radial range. The fit is satisfactory not only in the outer parts, where 
the density falls by nine orders of magnitude with respect to the center, but also in the inner 
regions. Depending on the initial virial parameter u, the end-products possess a density profile 
which, projected along the line-of-sight, may exhibit an R}'^ behavior with different values of 
n and yet it is well fitted by the f^^' models. In turii, this may be interpreted in the framework 
of the proposed weak homology of elliptical galaxies'^'. 

To some extent, the final anisotropy profiles for clumpy initial conditions are found to be 
sensitive to the detailed choice of initialization, in particular to the number and the size of the 
clumps. The agreement between end-products of the collapse and models (see Fig. ^ seems to 
be best for simulations with 10 to 20 clumps, for a clump size such that the sum of their volumes 
fills the sphere where their centers of mass lie. We should recall that both the low and the high 
number of clumps limits are those of an initial homogeneous condition, unfavorable to violent 
relaxation. 

At the level of phase space, we have performed two types of comparison, one involving the 





Figure 2: Comparison between the results of one simulation with u = 0.23, called C2_4 (starting from 20 cold 
clumps), and the best-fit /'"' model {v = 5/8, ^ — 5.4). Left set of four panels. Density as measured from the 
simulation (error bars) and model profile (top left). Residuals from the best-fit profile (top right). Anisotropy 
profile of the simulation (error bars) and profile (bottom left). Energy density distribution N{E) (bottom right). 
Right set of four panels. Final phase space density N{E, J^) (left column), compared with that of the best fitting 
Z*^"' model (right column). The model curve for N{E) and surface for N{E, J^) have been computed by a Monte 

Carlo sampling of phase space. 

energy density distribution N{E) and the other based on N{E, J^). The chosen normahzation 
factors are such that: M = f N{E)dE = J N{E,J^)dEdJ'^. In Fig. El we plot the final energy 
density distribution for a simulation run called C2_4 with respect to the predictions of the best- 
fit model identified from the study of the density and pressure anisotropy distributions. The 
agreement is good, especially for strongly bound particles. At the deeper level of N{E,J'^), 
simulations and models also agree quite well, as illustrated in the right set of four panels of 
Fig. 121 For the case shown, the distribution contour lines basically match in the range from 
Emin to E ^ —30 (in the relevant units); however, the theoretical model shows a peak located 
near the origin, not observed in the simulations. 
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